# *****************************************************************
# OVERVIEW ####
# *****************************************************************
# PSRM_REPLICATION.R
# Replication Code for:
# Are Rural Attitudes Just Republican?
# PSRM, 2023
# Authors: 
#   Jennifer Lin: jenniferlin2025@u.northwestern.edu
#   Kristin Lunz Trujillo: ktrujillo@hks.harvard.edu
# Created On: 2023 01 27

# ***The following script is intended to be an overview of all
#   the script files present in thsi directory and is intended for 
#   use with quick results. For more in-depth code explainations,
#   consult the respective code files.***

# *****************************************************************
# PACKAGES AND FUNCTIONS ####
# *****************************************************************

# Read in functions file for all the packages and 
#  functions for the analyses with the ANES Data

# *****************************************************************
# PACKAGES ####
# *****************************************************************

library(dplyr)       # For %>%, filter(), select(), mutate()
library(tidyr)       # For pivot_*() functions
library(haven)       # For read_dta()
library(srvyr)       # For as_survey()
library(survey)      # For survey summary statistics
library(broom)       # For tidy()
library(ggplot2)     # For plots
library(stargazer)   # For regression table exports
library(dotwhisker)  # For dot whisker plots in raw output
library(xtable)      # For table exports

# *****************************************************************
# FUNCTIONS ####
# *****************************************************************

# Below are a set of functions that are written by the authors 
#  that are used in the paper. This Script serves as the central
#. hub for locating all the custom functions that are not 
#  included in the packages previously loaded above.

## Cleaning Data ####

# Functions in this section address data cleaning. Since the 
#  process requires using many variables that have a similar 
#  coding structure, the functions help streamline the amount
#  of typing needed for this process.

# This is the general clean function to recode values that the 
#  ANES labels as missing (typically negative values) as NA. 
#  Cross-ref codebook before use.

clean_ANES <- function(var){
  dplyr::recode(
    var, 
    `-1` = NaN, `-2` = NaN, `-3` = NaN, `-4` = NaN,
    `-5` = NaN, `-6` = NaN, `-7` = NaN, `-8` = NaN,
    `-9` = NaN, `99` = NaN,
    .default = var)
} 

## Descriptive Statistics ####

# Functions in this section address the descriptive statistics
#  aspects of the paper.

# To generate survey means, the default for the 
#   `srvyr::survey_mean()` function has `na.rm = FALSE`. This
#   function resets the `na.rm` argument to be TRUE

svymNA <- function(var){
  survey_mean(var, na.rm = TRUE)
}

# In Supplemental Appendix C, we compute Chi-Squared statistics
#   for each of the attitudinal variables with the Party-Residence
#   variable as the predictor. To allow us to iterate across
#   all the variables in a seemless manner, the following function
#   computes survey weighted chi-squared statistics and prints 
#   the chi-squared statistic and its associated p-values.

chi_squ_svy <- function(data, columns, ...) {
  model <- lapply(as.list(columns), function(x) {
    svychisq(
      as.formula(paste0("~Rural_Party+", names(data)[x])), 
      data = data, statistic = c( "Chisq"), ...)
  })
  return(model)
}

## Regression Functions ####

# Since we have many outcome variables of interest in the main 
#   paper, it would be easier to iterate through them using a 
#   function as opposed to writing each model separately.
#   Therefore, the following functions accomplish the 
#   same purpose though using different models as 
#   described in the paper.

# Both functions allow us to iterate through all the models rather
#   effortlessly. The first function provides the model for
#   the main paper and tables in Supplemental Appendix D

issues_model <- function(data, columns, ...) {
  model <- lapply(as.list(columns), function(x) {
    svyglm(
      as.formula(
        paste0(
          names(data)[x],
          "~ Rural_Party + TRUMP + 
          ideo7 + FEMALE + income + educ + age + 
          rr_work + rr_slavery + rr_less + rr_harder +
          MINORITY + CHURCH")), 
      data = data, 
      ...)
  })
  return(model)
}

# The second function provides code for the models in Supplemetal 
#   Appendix E

issues_model_SI <- function(data, columns, ...) {
  model <- lapply(as.list(columns), function(x) {
    svyglm(
      as.formula(
        paste0(
          names(data)[x],
          "~ RURAL + PARTY + RURAL:PARTY + TRUMP + 
          ideo7 + FEMALE + income + educ + age + 
          rr_work + rr_slavery + rr_less + rr_harder +
          MINORITY + CHURCH")), 
      data = data, 
      ...)
  })
  return(model)
}

# In the Supplemental Appendix, we check whether respondents in the
#   cities differed substantially from suburban respondents.
# Here is the function for this model:
SC_model <- function(data, columns, ...) {
  model <- lapply(as.list(columns), function(x) {
    svyglm(
      as.formula(
        paste0(
          names(data)[x],
          "~ Place_Party + TRUMP + 
          ideo7 + FEMALE + income + educ + age +
          rr_work + rr_slavery + rr_less + rr_harder +
           MINORITY + CHURCH")), 
      data = data, 
      ...)
  })
  return(model)
}

# *****************************************************************
# LOAD DATA ####
# *****************************************************************

# All of the script files located in this replication directory 
#  requires that you have obtained a copy of the ANES 2020 data
#  from https://electionstudies.org/ and have loaded it into the
#  environment in your R session.

url <- "https://electionstudies.org/"
extension <- "wp-content/uploads/2022/02/"
filename <- "anes_timeseries_2020_stata_20220210.zip"
data_url <- paste0(url, extension, filename)

# Use BrowseURL to download the contents. Unzip them into a 
# folder called `data/ANES2020/` in the same file as this code.

browseURL(data_url)

# The below commands load the ANES data into your R Studio session
#  The ability to make use of this will depend on how you got the
#  original data files. To use these functions as is, simply place
#  the downloaded zip file into a folder called `data/` and name 
#  output of the zip as `ANES2020/`

data_file <- dir(
  here::here("data/ANES2020"),
  pattern    = ".+.dta$",    # Finds Stata Files
  full.names = TRUE,         # List Full File Path Names
  recursive  = TRUE)         # Repeat if needed


ANES_2020 <- read_dta(
  data_file
) %>% 
  zap_labels()

# *****************************************************************
# PREPARE DATA ####
# *****************************************************************

# The following code generates a cleaned version of the ANES Data
#   for use with the main paper. Key variables that we used are
#   cleaned and renamed in a more intuitive manner for ease of
#   reference.

ANES_Issues <- ANES_2020 %>% 
  mutate(
    # Clean and generate place of residence indicators
    residence       = clean_ANES(V202355),
    resid_Place     = case_when(
      residence == 1 ~ "Rural Area",
      residence == 2 ~ "Small Town",
      residence == 3 ~ "Suburb",
      residence == 4 ~ "City"
    ),
    resid_Place = factor(
      resid_Place,
      levels = c("City", "Suburb", "Small Town", "Rural Area")
    ),
    r_u_identity    = clean_ANES(V202356),
    resid_identity  = case_when(
      r_u_identity == 1 ~ "City Person",
      r_u_identity == 2 ~ "Suburb Person",
      r_u_identity == 3 ~ "Small Town Person",
      r_u_identity == 4 ~ "Country Person",
      r_u_identity == 5 ~ "Something else"
    ),
    resid_identity = factor(
      resid_identity,
      levels = c(
        "City Person", 
        "Suburb Person", 
        "Small Town Person", 
        "Country Person", 
        "Something Else")
    ),
    # Dichotomous variable to reflect Rural residence
    RURAL           = case_when(
      residence %in% c(1, 2) ~ TRUE,
      TRUE ~ FALSE
    ),
    # Generate Partisanship Variable
    pid7            = clean_ANES(V201231x),
    PARTY           = case_when(
      pid7 %in% c(1:3) ~ "Democrat",
      pid7  ==    4    ~ "Independent",
      pid7 %in% c(5:7) ~ "Republican",
      TRUE ~ NA_character_
    ),
    # Generate Place of Residence and Party Combination Variable
    Rural_Party = case_when(
      PARTY == "Democrat" & residence %in% c(1, 2) ~ 
        "Democrat-Rural",
      PARTY == "Democrat" & residence %in% c(3, 4) ~ 
        "Democrat-Urban",
      PARTY == "Republican" & residence %in% c(1, 2) ~ 
        "Republican-Rural",
      PARTY == "Republican" & residence %in% c(3, 4) ~ 
        "Republican-Urban"
    ),
    # Generate more nuanced place and party variable
    Place_Party = case_when(
      resid_Place == "City" & PARTY == "Democrat" ~ 
        "Democrat - City",
      resid_Place == "Suburb" & PARTY == "Democrat" ~ 
        "Democrat - Suburb",
      resid_Place == "Small Town" & PARTY == "Democrat" ~ 
        "Democrat - Small Town",
      resid_Place == "Rural Area" & PARTY == "Democrat" ~ 
        "Democrat - Rural Area",
      resid_Place == "City" & PARTY == "Republican" ~ 
        "Republican - City",
      resid_Place == "Suburb" & PARTY == "Republican" ~ 
        "Republican - Suburb",
      resid_Place == "Small Town" & PARTY == "Republican" ~ 
        "Republican - Small Town",
      resid_Place == "Rural Area" & PARTY == "Republican" ~ 
        "Republican - Rural Area"
    ),
    # Controls for Regression
    # Political Ideology
    ideo7           = clean_ANES(V201200),
    # Gender -- Indicator for Female
    gender          = clean_ANES(V201600),
    FEMALE          = case_when(
      gender == 2 ~ TRUE,
      TRUE ~ FALSE
    ),
    # Education
    educ            = clean_ANES(V201511x),
    # Race -- Indicator for Racial Minority
    race            = clean_ANES(V201549x),
    MINORITY        = case_when(
      race == 1 ~ FALSE,
      TRUE ~ TRUE
    ),
    # Income
    income          = clean_ANES(V201617x),
    # Age
    age             = clean_ANES(V201507x),
    # Church Attendance
    CHURCH          = clean_ANES(V201453),
    CHURCH          = 5 - CHURCH,
    # VOTE CHOICE in 2020 -- Indicators for Biden/Trump
    BIDEN = case_when(
      V202073 == 1 ~ TRUE,
      TRUE ~ FALSE
    ),
    TRUMP = case_when(
      V202073 == 2 ~ TRUE,
      TRUE ~ FALSE
    ),
    Pres_Vote = case_when(
      V202073 == 1 ~ "Joe Biden",
      V202073 == 2 ~ "Donald Trump"
    ),
    # Political Issue Attitudes
    # Allowing refugees who are fleeing war, persecution, or
    #   natural disasters in other countries to come to live 
    #   in the U.S.
    I_refugees = clean_ANES(V202236x),
    I_refugees = 8 - I_refugees,
    # Providing a path to citizenship for unauthorized 
    #   immigrants who obey the law, pay a fine, and
    #   pass security checks?
    I_citizenship = clean_ANES(V202242x),
    I_citizenship = 8 - I_citizenship,
    # Returning all unauthorized immigrants to their 
    #   native countries?
    I_deport = clean_ANES(V202245x),
    I_deport = 8 - I_deport,
    # Separating the children of detained immigrants, rather than 
    #   keeping them with their parents in adult detention centers?
    I_SPChild = clean_ANES(V202248x),
    I_SPChild = 8 - I_SPChild,
    # The government trying to reduce the difference in incomes 
    #   between the richest and poorest households?
    I_reduceIneq = clean_ANES(V202259x),
    I_reduceIneq = 8 - I_reduceIneq,
    # Affordable Care Act of 2010
    I_ACA = clean_ANES(V202328x),
    I_ACA = 8 - I_ACA,
    # Requiring children to be vaccinated in order to attend 
    #   public schools?
    I_reqVax = clean_ANES(V202331x),
    I_reqVax = 8 - I_reqVax,
    # Increased government regulation on businesses that 
    #   produce a great deal of greenhouse emissions linked 
    #   to climate change?
    I_RegGHG = clean_ANES(V202336x),
    I_RegGHG = 8 - I_RegGHG,
    # Requiring background checks for gun purchases at 
    #   gun shows or other private sales?
    I_bkgCheck = clean_ANES(V202341x),
    I_bkgCheck = 8 - I_bkgCheck,
    # Banning the sale of semi-automatic "assault-style" rifles?
    I_banAR = clean_ANES(V202344x),
    I_banAR = 8 - I_banAR,
    # A mandatory program where the government would buy back 
    #   semi-automatic assault-style rifles from citizens who 
    #   currently own them?
    I_buyBackAR = clean_ANES(V202347x),
    I_buyBackAR = 8 - I_buyBackAR,
    # Federal government should be doing more about 
    #   the opioid drug addiction issue
    I_Opioid = clean_ANES(V202350x),
    I_Opioid = 8 - I_Opioid,
    # The U.S. making free trade agreements with other countries?
    I_FreeTrade = clean_ANES(V202361x),
    I_FreeTrade = 8 - I_FreeTrade,
    # #stablishing a federal program that gives all citizens 
    #   $12,000 per year, provided they meet certain conditions?
    I_UBI = clean_ANES(V202376x),
    I_UBI = 8 - I_UBI,
    # Change in government spending to help people pay for health 
    #   insurance when they can’t pay for it all themselves?
    I_incHCspend = clean_ANES(V202380x),
    I_incHCspend = 8 - I_incHCspend,
    # Allowing transgender people to serve in the United States 
    #   Armed Forces?
    I_TransMil = clean_ANES(V202390x),
    I_TransMil = 8 - I_TransMil,
    # Requiring all people to show a government issued 
    #   photo ID when they vote?
    I_ReqIDVote = clean_ANES(V201359x),
    I_ReqIDVote = 8 - I_ReqIDVote,
    # Allowing convicted felons to vote once they complete 
    #   their sentence?
    I_FelonVote = clean_ANES(V201362x),
    I_FelonVote = 8 - I_FelonVote,
    # Elected officials restricting journalists' access to 
    #   information about government decision-making?
    I_JAccess = clean_ANES(V201375x),
    I_JAccess = 8 - I_JAccess,
    # Requiring employers to offer paid leave to parents of 
    #   new children?
    I_FamilyLeave = clean_ANES(V201405x),
    I_FamilyLeave = 8 - I_FamilyLeave,
    # Some people have proposed that the U.S. Constitution should 
    #   be changed so that the children of unauthorized immigrants 
    #   do not automatically get citizenship if they are born 
    #   in this country.
    I_birthright = clean_ANES(V201420x),
    I_birthright = 8 - I_birthright,
    # Building a wall on the U.S. border with Mexico?
    I_BorderWall = clean_ANES(V201426x),
    I_BorderWall = 8 - I_BorderWall,
    # Racial Resentment
    # Irish, Italian, Jewish and many other minorities overcame
    #   prejudice and worked their way up. Blacks should do 
    #   the same without any special favors.
    rr_work         = clean_ANES(V202300),
    # Generations of slavery and discrimination have created 
    #   conditions that make it difficult for blacks to work their 
    #   way out of the lower class.
    rr_slavery      = clean_ANES(V202301),
    # Over the past few years, blacks have gotten less 
    #   than they deserve.
    rr_less         = clean_ANES(V202302),
    # It's really a matter of some people not trying hard enough; 
    #   if blacks would only try harder they could be just as well 
    #   off as whites.
    rr_harder       = clean_ANES(V202303)
  )

# We utilize survey weights in the data analyses. Here, we use the
#   provided survey weight for the full sample pre-election
#   weight given that most of our variables originate from the
#   pre-election wave.

ANES_Survey <- ANES_Issues %>% 
  srvyr::as_survey(weight = V200010a)

# We can clear the workspace to just the objects that
#   we will use in the analyses.

rm(ANES_2020)

# *****************************************************************
# DESCRIPTIVE STATISTICS IN MAIN PAPER ####
# *****************************************************************

## Figure 1 -- Distribution of Respondents ####

# In Figure 1, we computed the number of respondents who belong
#   in each of the party-place of residence categories.

# To do this, we first generate a survey-weighted count of the
#   number of participants in each category
Rural_Party <- ANES_Survey %>% 
  group_by(Rural_Party) %>% 
  survey_count() %>% 
  filter(!is.na(Rural_Party))

# Then, we generate a plot to reflect the number of respondents
#   plus the survey-weighted standard errors generated from above
ggplot(
  Rural_Party, aes(x = Rural_Party, y = n, fill = Rural_Party))+
  geom_bar(
    stat = 'identity', 
    position = position_dodge(), color = "black")+
  geom_errorbar(
    aes(ymin = n-1.96*n_se, ymax = n+1.96*n_se), 
    width = 0.1, size = 1)+
  scale_fill_manual(
    name = "Residence + Party",
    values = c(
      "Democrat-Rural" = "#6baed6",
      "Democrat-Urban" = "#08519c",
      "Republican-Rural" = "#fb6a4a",
      "Republican-Urban" = "#cb181d"
    )
  )+
  xlab("Residence + Party Classification")+
  ylab("Weighted Number of Respondents")+
  labs(
    title = "Distribution of ANES Respondents",
    subtitle = "By Place of Residence and Party Classification",
    caption = 
      "Data: American National Elections Stidues (ANES) 2020"
  )+
  theme_bw()+
  theme(
    plot.title    = element_text(
      hjust = 0.5, size = 20, colour="black", face = "bold"),
    plot.subtitle = element_text(
      hjust = 0.5, size = 16, colour="black"),
    legend.title  = element_text(
      hjust = 0.5, size = 14, colour="black", face = "bold"),
    plot.caption  = element_text(
      size = 10, colour="black"),
    axis.title    = element_text(
      size = 14, colour="black", face = "bold"),
    axis.text     = element_text(
      size = 12, colour="black", angle = 0, hjust = 0.5),
    legend.position = "none"
  )


## Figure 2 ####

# In Figure 2, we compute the average ratings for participants
#   in each of the party-place of residence categories. 

# To do this, we 
#   - Compute by party-place categorization
#   - Generate a survey weighted mean to represent the average
#     response value for each group
#   - Compute a 95% confidence interval around the mean
issue_avg <- ANES_Survey %>% 
  group_by(Rural_Party) %>% 
  summarise_at(vars(starts_with("I_")), svymNA) %>% 
  filter(!is.na(Rural_Party)) %>% 
  reshape2::melt(by = c(Rural_Party), na.rm = TRUE) %>%
  mutate(valtype = ifelse(grepl("_se", variable), "se","mean")) %>%
  mutate(variable = gsub("_se", "", variable)) %>%
  pivot_wider(
    names_from = "valtype",
    values_from = "value") %>%
  mutate(
    lwr = mean - 1.96*se,
    upr = mean + 1.96*se,
    Issue_Desc = case_when(
      variable == "I_ACA" ~ "Approve Affordable Care Act",
      variable == "I_banAR" ~ "Ban Assault Style Rifles",
      variable == "I_birthright" ~ "End Birthright Citizenship",
      variable == "I_bkgCheck" ~ 
        "Background Check for\n Gun Purchases",
      variable == "I_BorderWall" ~ "Build Wall on Southern Border",
      variable == "I_buyBackAR" ~ "Assault Rifle Buy Back",
      variable == "I_citizenship" ~ 
        "Provide a Path to Citizenship",
      variable == "I_deport" ~ 
        "Deporting Immigrants\n to Native Country",
      variable == "I_FamilyLeave" ~ "Provide Paid Family Leave",
      variable == "I_FelonVote" ~ "Allowing Felons to Vote",
      variable == "I_FreeTrade" ~ "Allowing Free Trade Agreements",
      variable == "I_incHCspend" ~ 
        "Increasing Spending\n to Cover Health Care",
      variable == "I_JAccess" ~ "Transparency for Journalists",
      variable == "I_Opioid" ~ 'Address Opioid Epidemic',
      variable == "I_reduceIneq" ~ "Reduce Income Inequality",
      variable == "I_refugees" ~ "Allow Refugees to Come to US",
      variable == "I_RegGHG" ~ 
        "Regulate Greenhouse\n Gas Emissions",
      variable == "I_ReqIDVote" ~ "Require ID to Vote",
      variable == "I_reqVax" ~ "Require Vaccinations",
      variable == "I_SPChild" ~ 
        "Separating Children from\n Parents at Border",
      variable == "I_TransMil" ~ 
        "Allow Transgender People\n to Serve in Military",
      variable == "I_UBI" ~ "Provide Citizens 12K A Year"
    )
  )

# We can generate the plot using the information above.
ggplot(
  issue_avg, 
  aes(x = Issue_Desc, y = mean, color = Rural_Party))+
  geom_pointrange(
    aes(ymin = lwr, ymax = upr),
    position = position_dodge(.7))+
  scale_color_manual(
    name = "Residence + Party",
    values = c(
      "Democrat-Rural" = "#6baed6",
      "Democrat-Urban" = "#08519c",
      "Republican-Rural" = "#fb6a4a",
      "Republican-Urban" = "#cb181d"
    )
  )+
  xlab("Issues")+
  ylab("Mean")+
  labs(
    title = "Posiions on Core Issues",
    subtitle = "By Party and Place of Residence",
    caption = "Data: American National Elections Studies 2020"
  )+
  coord_flip()+
  guides(color = guide_legend(ncol = 2))+
  theme_bw()+
  theme(
    plot.title    = element_text(
      hjust = 0.5, size = 20, colour="black", face = "bold"),
    plot.subtitle = element_text(
      hjust = 0.5, size = 16, colour="black"),
    legend.title  = element_text(
      hjust = 0.5, size = 14, colour="black", face = "bold"),
    plot.caption  = element_text(
      size = 10, colour="black"),
    axis.title    = element_text(
      size = 14, colour="black", face = "bold"),
    axis.text     = element_text(
      size = 12, colour="black", angle = 0, hjust = 0.5),
    legend.position = "bottom"
  )

# For ease of use in the Chi-Squared tests below, we pull a vector
#   representing the descriptions for the issues for use later on
issues <- unique(issue_avg$Issue_Desc)

# *****************************************************************
# CHI-SQUARED TESTS IN SUPPLEMENTAL APPENDIX ####
# *****************************************************************

## Chi-Squared tests in Supplemental Appendix C ####

# To generate a Chi-Squared test for the issue positions by our
#   party-place of residence variable, we 
#   - Generate a second survey object so that the variables for the
#     issue positions are treated categorically rather than
#     continuously. 
#   - From here, we can compute the Chi-Squared statistic for each
#     issue position.
ANES_chi <- ANES_Issues %>% 
  filter(!is.na(Rural_Party)) %>% 
  mutate_at(vars(starts_with("I_")), as.character) %>% 
  srvyr::as_survey(weight = V200010a)

# We extract a data frame that consists of the variables that 
#   measure issue attitudes, analyzed above.
issue_vars <- ANES_Issues %>% 
  select(starts_with("I_")) 

# We compute a chi-squared statistic for each of the variables
#   by iterating through the variables for effeciency.
chi_issues <- chi_squ_svy(issue_vars, 1:22, design = ANES_chi)

# Create empty vector to store Chi-Squared estimate
x_sq <- c()
# Create empty vector to store p-value
x_p <- c()

# Use the following for loop to parse the chi-squared statistic
#   and p value from each of the analyses stored in the resulting
#   list from above.
for (i in 1:length(chi_issues)) {
  x_sq[i] <- chi_issues[[i]]$statistic
  x_p[i] <- chi_issues[[i]]$p.value
}

# From here, we can export the data as one concice table, as seen
#   in Supplemental Appendix C.
chi_sq_output <- data.frame(
  Issues = issues, 
  `X-Squared` = x_sq,
  `p.value` = x_p
) %>% 
  mutate(
    signif = case_when(
      p.value <= 0.001 ~ "***",
      p.value > 0.001 & p.value <= 0.05 ~ "**",
      p.value > 0.05 & p.value <= 0.1 ~ "*",
      TRUE ~ ""
    ),
    `Residence and Party` = paste0(
      round(`X.Squared`, 2), signif
    )
  ) %>% 
  select(Issues, `Residence and Party`)

# *****************************************************************
# INFERENTIAL STATISTICS IN MAIN PAPER ####
# *****************************************************************

# The following code generates the tables and figures that are 
#   present in the main paper.

# To start, we extract the variables that we want to regress. 
vars_to_regress <- ANES_Issues %>%
  select(starts_with("I_"))

# To run the regression using survey weights, we declare the data
#   as a survey object and remove Independents
ANES_regress_Data <- ANES_Issues %>% 
  filter(PARTY != "Independent") %>% 
  srvyr::as_survey(weight = V200010a)

# We iterate through each of the variables and apply the model 
#   specified in the paper and in the `issues_model()` function.
I_results <- issues_model(
  vars_to_regress, 1:22, design = ANES_regress_Data)

# Create a vector of terms to label predictors in the models to 
#   use with the regression table exports.
terms <- c(
  "Democrat-Urban",
  "Republican-Rural",
  "Republican-Urban",
  "Trump Supporter",
  "Ideology",
  "Female",
  "Income",
  "Education",
  "Age",
  "Racial Resentment -- Work without Favors",
  "Racial Resentment -- History of Slavery",
  "Racial Resentment -- Blacks get less",
  "Racial Resentment -- Try Harder",
  "Racial Minority",
  "Church"
)

# We can plot the results in a graph, as shown in Figure 3 in the 
#   main paper

# To start, create an empty list to store results from each of the
#   regression models
plot_df <- list()

# To get the estimates and standard erros to plot, we pull these#
#   details from our list of regression models, iterating
#   trhough each model one by one
for (i in 1:length(I_results)) {
  plot_df[[i]] <- tidy(I_results[[i]])
  plot_df[[i]]$variable <- names(vars_to_regress[i])
}

# We then clean up the resulting data frame, renaming the 
#   varaible names as needed
plot_df_clean <- bind_rows(plot_df) %>% 
  filter(grepl("Rural_Party", term)) %>% 
  mutate(
    Issue_Desc = case_when(
      variable == "I_refugees" ~ "Allow Refugees to Come to US",
      variable == "I_citizenship" ~ "Provide a Path to Citizenship",
      variable == "I_deport" ~ "Deporting Immigrants to Native Country",
      variable == "I_SPChild" ~ "Separating Children from Parents at Border",
      variable == "I_reduceIneq" ~ "Reduce Income Inequality",
      variable == "I_ACA" ~ "Approve Affordable Care Act",
      variable == "I_reqVax" ~ "Require Vaccinations",
      variable == "I_RegGHG" ~ "Regulate Greenhouse Gas Emissions",
      variable == "I_bkgCheck" ~ "Background Check for Gun Purchases",
      variable == "I_banAR" ~ "Ban Assault Style Rifles",
      variable == "I_buyBackAR" ~ "Assault Rifle Buy Back",
      variable == "I_Opioid" ~ 'Address Opioid Epidemic',
      variable == "I_FreeTrade" ~ "Allowing Free Trade Agreements",
      variable == "I_UBI" ~ "Provide Citizens 12K A Year",
      variable == "I_incHCspend" ~ "Increasing Spending to Cover Health Care",
      variable == "I_TransMil" ~ "Allow Transgender People to Serve in Military",
      variable == "I_ReqIDVote" ~ "Require ID to Vote",
      variable == "I_FelonVote" ~ "Allowing Felons to Vote",
      variable == "I_JAccess" ~ "Transparency for Journalists",
      variable == "I_FamilyLeave" ~ "Provide Paid Family Leave",
      variable == "I_birthright" ~ "End Birthright Citizenship",
      variable == "I_BorderWall" ~ "Build Wall on Southern Border"
    ),
    Issue_Desc = factor(Issue_Desc, levels = unique(Issue_Desc)),
    term = case_when(
      term == "Rural_PartyDemocrat-Urban" ~ "Democrat-Urban",
      term == "Rural_PartyRepublican-Rural" ~ "Republican-Rural",
      term == "Rural_PartyRepublican-Urban" ~ "Republican-Urban"
    ),
    term = factor(
      term,
      levels = c(
        "Democrat-Urban", "Republican-Rural", "Republican-Urban"
      )
    )
  )

# Finally, we generate the plot
ggplot(plot_df_clean, aes(x = term, y = estimate))+
  geom_pointrange(
    aes(
      ymin = estimate-1.96*std.error, 
      ymax = estimate+1.96*std.error)
  )+
  geom_hline(
    yintercept = 0, color = "grey75", linetype = "dashed")+
  facet_wrap(~Issue_Desc, ncol = 2)+
  coord_flip()+
  labs(
    title = "Favor or Oppose?",
    subtitle = 
      "Evaluations of Core Issues by Place of Residence and Party",
    caption = 
      "Data: American National Elections Studies (ANES) 2020"
  )+
  ylab("Estimate")+
  xlab("Key Predictors")+
  theme_bw()+
  theme(
    plot.title    = element_text(
      hjust = 0.5, size = 20, colour="black", face = "bold"),
    plot.subtitle = element_text(
      hjust = 0.5, size = 16, colour="black"),
    legend.title  = element_text(
      hjust = 0.5, size = 14, colour="black", face = "bold"),
    plot.caption  = element_text(size = 10, colour="black"),
    axis.title    = element_text(size = 14, colour="black"),
    axis.text.x   = element_text(
      size = 12, colour="black", angle = 0, hjust = 0.5)
  )  

# *****************************************************************
# INFERENTIAL STATISTICS IN SUPPLEMENTAL APPENDIX ####
# *****************************************************************

# The following code generates the tables and figures that are 
#   present in the Supplemental Appendix.

# To start, we extract the variables that we want to regress. As
#   mentioned in the paper, the selection of variables is based on
#   significant differences within each party by urban-rural 
#   factors. This leads to the following variables:
vars_to_regress <- ANES_Issues %>%
  select(starts_with("I_"))

# To run the regression using survey weights, we declare the data
#   as a survey object and remove Independents
ANES_regress_Data <- ANES_Issues %>%
  filter(PARTY != "Independent") %>%
  srvyr::as_survey(weight = V200010a)

# We iterate through each of the variables and apply the model 
#   specified in the paper and in the `issues_model_si()` 
#   function.
I_results <- issues_model_SI(
  vars_to_regress, 1:22, design = ANES_regress_Data)

# Create a vector of terms to label predictors in the models to 
#   use with the regression table exports.
terms <- c(
  "Rural Resident",
  "Republican",
  "Trump Supporter",
  "Ideology",
  "Female",
  "Income",
  "Education",
  "Age",
  "Racial Resentment -- Work without Favors",
  "Racial Resentment -- History of Slavery",
  "Racial Resentment -- Blacks get less",
  "Racial Resentment -- Try Harder",
  "Racial Minority",
  "Church",
  "Rural Resident X Republican"
)

# We can plot the results in a graph, as shown in Supplemental
#   Appendix E.

# To start, create an empty list to store results from each of the
#   regression models
plot_df <- list()

# To get the estimates and standard erros to plot, we pull these#
#   details from our list of regression models, iterating
#   trhough each model one by one
for (i in 1:length(I_results)) {
  plot_df[[i]] <- tidy(I_results[[i]])
  plot_df[[i]]$variable <- names(vars_to_regress[i])
}

# We then clean up the resulting data frame, renaming the 
#   varaible names as needed
plot_df_clean <- bind_rows(plot_df) %>% 
  filter(term %in% c(
    "RURALTRUE", "PARTYRepublican", "RURALTRUE:PARTYRepublican")) %>% 
  mutate(
    Issue_Desc = case_when(
      variable == "I_refugees" ~ "Allow Refugees to Come to US",
      variable == "I_citizenship" ~ "Provide a Path to Citizenship",
      variable == "I_deport" ~ "Deporting Immigrants to Native Country",
      variable == "I_SPChild" ~ "Separating Children from Parents at Border",
      variable == "I_reduceIneq" ~ "Reduce Income Inequality",
      variable == "I_ACA" ~ "Approve Affordable Care Act",
      variable == "I_reqVax" ~ "Require Vaccinations",
      variable == "I_RegGHG" ~ "Regulate Greenhouse Gas Emissions",
      variable == "I_bkgCheck" ~ "Background Check for Gun Purchases",
      variable == "I_banAR" ~ "Ban Assault Style Rifles",
      variable == "I_buyBackAR" ~ "Assault Rifle Buy Back",
      variable == "I_Opioid" ~ 'Address Opioid Epidemic',
      variable == "I_FreeTrade" ~ "Allowing Free Trade Agreements",
      variable == "I_UBI" ~ "Provide Citizens 12K A Year",
      variable == "I_incHCspend" ~ "Increasing Spending to Cover Health Care",
      variable == "I_TransMil" ~ "Allow Transgender People to Serve in Military",
      variable == "I_ReqIDVote" ~ "Require ID to Vote",
      variable == "I_FelonVote" ~ "Allowing Felons to Vote",
      variable == "I_JAccess" ~ "Transparency for Journalists",
      variable == "I_FamilyLeave" ~ "Provide Paid Family Leave",
      variable == "I_birthright" ~ "End Birthright Citizenship",
      variable == "I_BorderWall" ~ "Build Wall on Southern Border"
    ),
    Issue_Desc = factor(Issue_Desc, levels = unique(Issue_Desc)),
    term = case_when(
      term == "RURALTRUE" ~ "Rural Resident",
      term == "PARTYRepublican" ~ "Republican",
      term == "RURALTRUE:PARTYRepublican" ~ "Rural X Republican"
    ),
    term = factor(
      term,
      levels = c(
        "Rural Resident", "Republican", "Rural X Republican"
      )
    )
  )

# Finally, we generate the plot
ggplot(plot_df_clean, aes(x = term, y = estimate))+
  geom_pointrange(
    aes(
      ymin = estimate-1.96*std.error, 
      ymax = estimate+1.96*std.error)
  )+
  geom_hline(
    yintercept = 0, 
    color = "grey75", linetype = "dashed")+
  facet_wrap(~Issue_Desc, ncol = 2)+
  coord_flip()+
  labs(
    title = "Favor or Oppose?",
    subtitle = 
      "Evaluations of Core Issues by Place of Residence and Party",
    caption = 
      "Data: American National Elections Studies (ANES) 2020"
  )+
  ylab("Estimate")+
  xlab("Key Predictors")+
  theme_bw()+
  theme(
    plot.title    = element_text
    (hjust = 0.5, size = 20, colour="black", face = "bold"),
    plot.subtitle = element_text(
      hjust = 0.5, size = 16, colour="black"),
    legend.title  = element_text
    (hjust = 0.5, size = 14, colour="black", face = "bold"),
    plot.caption  = element_text(
      size = 10, colour="black"),
    axis.title    = element_text(
      size = 14, colour="black"),
    axis.text.x   = element_text(
      size = 12, colour="black", angle = 0, hjust = 0.5)
  )  

# *****************************************************************
# CHECK SUBURB AND CITY DIFFERENCES ####
# *****************************************************************

# We can run a regression model using the model predefined
#   in the Functions.R file called "SC_model".
I_results <- SC_model(
  vars_to_regress, 1:22, design = ANES_regress_Data)

# Then, we can generate a plot that appears in the page.

# To start, create an empty list to store results from each of the
#   regression models
plot_df <- list()

# To get the estimates and standard erros to plot, we pull these#
#   details from our list of regression models, iterating
#   trhough each model one by one
for (i in 1:length(I_results)) {
  plot_df[[i]] <- tidy(I_results[[i]])
  plot_df[[i]]$variable <- names(vars_to_regress[i])
}

# We then clean up the resulting data frame, renaming the 
#   varaible names as needed
plot_df_clean <- bind_rows(plot_df) %>% 
  filter(grepl("Place_Party", term)) %>% 
  mutate(
    Issue_Desc = case_when(
      variable == "I_refugees" ~ "Allow Refugees to Come to US",
      variable == "I_citizenship" ~ "Provide a Path to Citizenship",
      variable == "I_deport" ~ "Deporting Immigrants to Native Country",
      variable == "I_SPChild" ~ "Separating Children from Parents at Border",
      variable == "I_reduceIneq" ~ "Reduce Income Inequality",
      variable == "I_ACA" ~ "Approve Affordable Care Act",
      variable == "I_reqVax" ~ "Require Vaccinations",
      variable == "I_RegGHG" ~ "Regulate Greenhouse Gas Emissions",
      variable == "I_bkgCheck" ~ "Background Check for Gun Purchases",
      variable == "I_banAR" ~ "Ban Assault Style Rifles",
      variable == "I_buyBackAR" ~ "Assault Rifle Buy Back",
      variable == "I_Opioid" ~ 'Address Opioid Epidemic',
      variable == "I_FreeTrade" ~ "Allowing Free Trade Agreements",
      variable == "I_UBI" ~ "Provide Citizens 12K A Year",
      variable == "I_incHCspend" ~ "Increasing Spending to Cover Health Care",
      variable == "I_TransMil" ~ "Allow Transgender People to Serve in Military",
      variable == "I_ReqIDVote" ~ "Require ID to Vote",
      variable == "I_FelonVote" ~ "Allowing Felons to Vote",
      variable == "I_JAccess" ~ "Transparency for Journalists",
      variable == "I_FamilyLeave" ~ "Provide Paid Family Leave",
      variable == "I_birthright" ~ "End Birthright Citizenship",
      variable == "I_BorderWall" ~ "Build Wall on Southern Border"
    ),
    Issue_Desc = factor(Issue_Desc, levels = unique(Issue_Desc)),
    term = case_when(
      term == "Place_PartyDemocrat - Suburb" ~ "Democrat - Suburb",
      term == "Place_PartyDemocrat - Small Town" ~ "Democrat - Small Town",
      term == "Place_PartyDemocrat - Rural Area" ~ "Democrat - Rural Area",
      term == "Place_PartyRepublican - City" ~ "Republican - City",
      term == "Place_PartyRepublican - Suburb" ~ "Republican - Suburb",
      term == "Place_PartyRepublican - Small Town" ~ "Republican - Small Town",
      term == "Place_PartyRepublican - Rural Area" ~ "Republican - Rural Area"
    ),
    term = factor(
      term,
      levels = c(
        "Democrat - Suburb", "Democrat - Small Town", "Democrat - Rural Area",
        "Republican - City", "Republican - Suburb", "Republican - Small Town",
        "Republican - Rural Area"
      )
    )
  ) %>% 
  filter(grepl("Suburb|City", term))

# Finally, we generate the plot
ggplot(plot_df_clean, aes(x = term, y = estimate))+
  geom_pointrange(
    aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error)
  )+
  geom_hline(yintercept = 0, color = "grey75", linetype = "dashed")+
  facet_wrap(~Issue_Desc, ncol = 2)+
  coord_flip()+
  labs(
    title = "Favor or Oppose?",
    subtitle = "Evaluations of Core Issues by Place of Residence and Party",
    caption = "Data: American National Elections Studies (ANES) 2020"
  )+
  ylab("Estimate")+
  xlab("Key Predictors")+
  theme_bw()+
  theme(
    plot.title    = element_text(hjust = 0.5, size = 20, colour="black", face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 16, colour="black"),
    legend.title  = element_text(hjust = 0.5, size = 14, colour="black", face = "bold"),
    plot.caption  = element_text(size = 10, colour="black"),
    axis.title    = element_text(size = 14, colour="black"),
    axis.text.x   = element_text(size = 12, colour="black", angle = 0, hjust = 0.5)
  ) 
